MODULE pfmm_main

   USE fmm_global_consts
   USE fmm_global_types
   USE pfmm_global_types

   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: pfmm_initial,                      &
             pfmm_final,                        &
             pfmm_get_J_matrix

   ! The multipole potential
   REAL(REALK), POINTER, SAVE :: Vff(:,:)
   ! The (modified) multipole expansions
   TYPE(raw_mm_data), SAVE :: LHS_mms, RHS_mms
   TYPE(T_vector_map), POINTER, SAVE :: Tvec_map(:)

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_initial(LMAX)

      USE fmm_scheme_module, ONLY: fmm_init_scheme
      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: LMAX

      CALL fmm_init_scheme(PBC_FMM,LMAX)

   END SUBROUTINE pfmm_initial

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_final
      IMPLICIT NONE
      CALL pfmm_free_memory
   END SUBROUTINE pfmm_final

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_init_qlm_and_V(scheme,dens)

      USE pfmm_qlm_builder,    ONLY: pfmm_get_raw_qlm
      USE fmm_aux_qlm_builder, ONLY: fmm_get_aux_qlm

      IMPLICIT NONE
      TYPE(fmm_scheme),   INTENT(IN) :: scheme
      REAL(REALK),        INTENT(IN) :: dens(:,:,:)

      NULLIFY (Vff)
      NULLIFY (Tvec_map)
      ! Get basic multipole data
      CALL pfmm_get_raw_qlm(scheme,dens,LHS_mms,RHS_mms,Tvec_map)

      ! Now ensure appropriate prefactors, normalisation and packing
      CALL fmm_get_aux_qlm(scheme,LHS_mms,RHS_mms)
      ! Allocate the far field potential
      CALL pfmm_allocate_Vff(scheme)

   END SUBROUTINE pfmm_init_qlm_and_V

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_free_memory

      USE fmm_qlm_builder, ONLY: fmm_deallocate_qlm
      IMPLICIT NONE
      IF (ASSOCIATED(Vff)) DEALLOCATE(Vff)
      NULLIFY (Vff)
      CALL fmm_deallocate_qlm(LHS_mms,RHS_mms)
      IF (ASSOCIATED(Tvec_map)) DEALLOCATE(Tvec_map)
      NULLIFY (Tvec_map)

   END SUBROUTINE pfmm_free_memory

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_allocate_Vff(scheme)

      IMPLICIT NONE
      TYPE(fmm_scheme), INTENT(IN) :: scheme
      INTEGER(INTK) :: lm_dim, mms_dim, alloc_error

      IF (.NOT.ASSOCIATED(LHS_mms%paras)) STOP 'mms ptrs not set in fmm_main!'
      IF (ASSOCIATED(Vff)) CALL fmm_quit('Vff should NOT be allocated already!')

      mms_dim = SIZE(LHS_mms%paras)
      ! Note we shouldn't use this array for translated, BOXED potentials
      lm_dim = (1+ scheme%raw_LMAX)**2

      WRITE(IO1,*) 'Vff: Attempting to allocate',  &
                     MAX(1,lm_dim*mms_dim*8/1000000), 'MB of memory...'
      ALLOCATE (Vff(lm_dim,mms_dim), STAT=alloc_error)
      IF (alloc_error /= 0) WRITE(IO1,*) '... Failed!'

      ! Must zero out since Vff is built additively in general
      Vff(:,:) = zero

   END SUBROUTINE pfmm_allocate_Vff

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_get_J_matrix(dens,lattice_consts,J_matrix)

      USE fmm_scheme_module, ONLY: fmm_get_scheme
      USE pfmm_Vff_driver,   ONLY: pfmm_get_Vff
      USE pfmm_J_builder,    ONLY: pfmm_get_J_from_V

      IMPLICIT NONE
      REAL(REALK), INTENT(IN)  :: dens(:,:,:)
      REAL(REALK), INTENT(IN)  :: lattice_consts(3)
      REAL(REALK), INTENT(OUT) :: J_matrix(:,:,:)

      TYPE(fmm_scheme), POINTER :: scheme
      REAL(REALK) :: energy, T0, T1
integer(intk) :: i

      CALL CPU_TIME(T0)
      CALL fmm_get_scheme(scheme)

      ! We only have the density on RHS when getting J-matrix
      scheme%LHS_dens = .FALSE.
      scheme%RHS_dens = .TRUE.
      scheme%pack_LHS = .FALSE.
      scheme%pack_RHS = .TRUE.

      ! Prepare moments and allocate potential
      CALL pfmm_init_qlm_and_V(scheme,dens)
print *, "LHS qlm_T:"
do i = 1, size(LHS_mms%paras)
    print *, i, LHS_mms%qlm_T(1,i)
end do
print *, "RHS qlm_T:"
do i = 1, size(RHS_mms%paras)
    print *, i, RHS_mms%qlm_T(1,i)
end do

      ! Get potential
      CALL pfmm_get_Vff(scheme,lattice_consts,LHS_mms%paras,RHS_mms,Vff)

      ! Get J-matrix
      CALL pfmm_get_J_from_V(scheme,Vff,Tvec_map,LHS_mms,J_matrix,energy)
      WRITE(IO1,'(X,A," = ",E20.12)') "Total Coulomb energy =", energy
print *, "Total Coulomb energy =", energy

      CALL pfmm_free_memory
      CALL CPU_TIME(T1)
      CALL fmm_TIMTXT('>>> TIME USED in fmm_get_J_matrix', T1-T0, IO1)

   END SUBROUTINE pfmm_get_J_matrix

!-------------------------------------------------------------------------------

END MODULE pfmm_main

